'''
Created on Aug 2, 2012

@author: Xindong
'''


from shortcuts import execute_sql, execute_select, log_info
import math


#    +---------+---------+ lat_max
#    |         |         |
#    |    C    |    D    |
#    |         |         |
#    +---------+---------+ lat_mid
#    |         |         |
#    |    A    |    B    |
#    |         |         |
#    +---------+---------+ lat_min
# lon_min   lon_mid   lon_max
# 
# A = [lat_min, lat_mid) x [lon_min, lon_mid)
# B = [lat_min, lat_mid) x [lon_mid, lon_max)
# C = [lat_mid, lat_max) x [lon_min, lon_mid)
# D = [lat_mid, lat_max) x [lon_mid, lon_max)
# for the entire usa
#TILE_LAT_MIN = +2500000
#TILE_LAT_MAX = +5000000
#TILE_LON_MIN = -12500000
#TILE_LON_MAX = -6500000


''' For the moment, we dont create a real map database to record node distribution
    but the put tile_id in the node table and use like statement to  
'''

TILE_LAT_MIN = +3930000 #+3936329
TILE_LAT_MAX = +4190000 #+4182157
TILE_LON_MIN = -7520000 #-7518232
TILE_LON_MAX = -7190000 #-7195388

##london
#TILE_LAT_MIN = +5100000
#TILE_LAT_MAX = +5200000
#TILE_LON_MIN = -62000
#TILE_LON_MAX = +26000
#
##sweden
#TILE_LAT_MIN = +5200000
#TILE_LAT_MAX = +7000000
#TILE_LON_MIN = +1000000
#TILE_LON_MAX = +2500000


MIN_LAT_LON_SCALE = 50     # 60 m
TILE_SCALE = max( TILE_LAT_MAX - TILE_LAT_MIN, TILE_LON_MAX - TILE_LON_MIN ) / MIN_LAT_LON_SCALE
#MAX_TILE_LEVEL = LOG2(TILE_SCALE)
MAX_TILE_LEVEL = TILE_SCALE.bit_length() - 1


''' We are so happy to find that 0.00001 degree is about 1.1m in most part of the world
        delta_lat 0.00001 == delta_dis 1.1m
        delta_lon 0.00001 == delta_dis 1.1m * (1-lat*lat/8100)
    So we will keep the precision 0.00001 for lat/lon values
    Multiply by 100000, we have:
        delta_lat 1 == delta_dis 1.1m
        delta_lon 1 == delta_dis 1.1m * (1-lat*lat/8100/1e+10)
    Inverse:
        delta_dis 100 m == delta_lat 91
        delta_dis 100 m == delta_lon 91 / (1-lat*lat/8100/1e+10)
    Precision in most of the case
        with an error < 1.0m for short distances
        and an error < 6.0% for long distances
'''


def tiles_covering(rect, tile_rect = None, level = None):
    if tile_rect == None:
        tile_rect = TILE_LAT_MIN, TILE_LAT_MAX, TILE_LON_MIN, TILE_LON_MAX
        level = 0
    
    tile_lat_min, tile_lat_max, tile_lon_min, tile_lon_max = tile_rect
    lat_min, lat_max, lon_min, lon_max = rect
    lat_min, lat_max = max(lat_min, tile_lat_min), min(lat_max, tile_lat_max)
    lon_min, lon_max = max(lon_min, tile_lon_min), min(lon_max, tile_lon_max)
    
    if lat_min >= lat_max or lon_min >= lon_max:
        return []
    
    if tile_rect == rect or level == MAX_TILE_LEVEL:
        return ["%"]
    
    tile_list = []
    
    tile_lat_mid = (tile_lat_min + tile_lat_max) >> 1
    tile_lon_mid = (tile_lon_min + tile_lon_max) >> 1
    
    tile_rect = tile_lat_min, tile_lat_mid, tile_lon_min, tile_lon_mid
    rect = min(lat_min, tile_lat_mid), min(lat_max, tile_lat_mid), min(lon_min, tile_lon_mid), min(lon_max, tile_lon_mid)
    sub_tile_list = tiles_covering(rect, tile_rect, level + 1)
    for tile_id in sub_tile_list:
        tile_list.append("A" + tile_id)
    
    tile_rect = tile_lat_min, tile_lat_mid, tile_lon_mid, tile_lon_max
    rect = min(lat_min, tile_lat_mid), min(lat_max, tile_lat_mid), max(lon_min, tile_lon_mid), max(lon_max, tile_lon_mid)
    sub_tile_list = tiles_covering(rect, tile_rect, level + 1)
    for tile_id in sub_tile_list:
        tile_list.append("B" + tile_id)
        
    tile_rect = tile_lat_mid, tile_lat_max, tile_lon_min, tile_lon_mid
    rect = max(lat_min, tile_lat_mid), max(lat_max, tile_lat_mid), min(lon_min, tile_lon_mid), min(lon_max, tile_lon_mid)
    sub_tile_list = tiles_covering(rect, tile_rect, level + 1)
    for tile_id in sub_tile_list:
        tile_list.append("C" + tile_id)
   
    tile_rect = tile_lat_mid, tile_lat_max, tile_lon_mid, tile_lon_max
    rect = max(lat_min, tile_lat_mid), max(lat_max, tile_lat_mid), max(lon_min, tile_lon_mid), max(lon_max, tile_lon_mid)
    sub_tile_list = tiles_covering(rect, tile_rect, level + 1)
    for tile_id in sub_tile_list:
        tile_list.append("D" + tile_id)
    
    return tile_list


def calculate_tile_id(latitude, longitude):
    tile_lat_min = TILE_LAT_MIN
    tile_lat_max = TILE_LAT_MAX
    tile_lon_min = TILE_LON_MIN
    tile_lon_max = TILE_LON_MAX
    leaf_id = ""
    for _i in range(MAX_TILE_LEVEL):
        tile_lat_mid = (tile_lat_min + tile_lat_max) >> 1
        tile_lon_mid = (tile_lon_min + tile_lon_max) >> 1
        if latitude < tile_lat_mid and longitude < tile_lon_mid:
            tile_lat_max, tile_lon_max = tile_lat_mid, tile_lon_mid
            leaf_id += "A"
        elif latitude < tile_lat_mid and longitude >= tile_lon_mid:
            tile_lat_max, tile_lon_min = tile_lat_mid, tile_lon_mid
            leaf_id += "B"
        elif latitude >= tile_lat_mid and longitude < tile_lon_mid:
            tile_lat_min, tile_lon_max = tile_lat_mid, tile_lon_mid
            leaf_id += "C"
        elif latitude >= tile_lat_mid and longitude >= tile_lon_mid:
            tile_lat_min, tile_lon_min = tile_lat_mid, tile_lon_mid
            leaf_id += "D"
    return leaf_id


def update_node_tiles(dataset_name):
    log_info("tiles calculation for " + dataset_name)
    nodes = execute_select("""
        SELECT dataset_id, node_id, latitude, longitude
        FROM %s_graph_nodes
    """ % dataset_name)
    for dataset_id, node_id, latitude, longitude in nodes:
        tile_id = calculate_tile_id(latitude, longitude)
        execute_sql("""
            UPDATE %s_graph_nodes SET tile_id = '%s' WHERE dataset_id = %d AND node_id = %d
        """ % (dataset_name, tile_id, dataset_id, node_id))
    log_info("tiles calculation end")


def find_neighbours(latitude, longitude, radius, from_where = "gtfs"):
    ''' @var radius: in meter
        neighbour area:
                  #
              # # # # #
            # # # # # # #
            # # # # # # #
          # # # # o # # # #
            # # # # # # #
            # # # # # # #
              # # # # #
                  #
        Manhattan distance < 5/4
    '''
    rect = circle2tile(latitude, longitude, radius)
    tile_list = tiles_covering(rect)
    neighbours = []
    if len(tile_list) > 0: 
        like_list = " or ".join([("tile_id like '%s'" % tile_id) for tile_id in tile_list])
        results = execute_select("""
            SELECT dataset_id, node_id, latitude, longitude
            FROM %s_graph_nodes
            WHERE %s
        """ % (from_where, like_list))
        for dataset_id, node_id, node_lat, node_lon in results:
            if abs(latitude - int(node_lat)) + abs(longitude - int(node_lon)) <= 5 * radius / 4:
                neighbours.append((int(dataset_id), int(node_id)))
    return neighbours


##################################################################################################


def circle2tile(latitude, longitude, radius):
    ''' delta_lat = int(radius / 1.1)
        delta_lon = int(radius / 1.1 / (1-latitude*latitude/1e+10/8100))
    '''
    delta_lat = radius * 10 / 11
    delta_lon = delta_lat * 81000000000000 / (81000000000000 - latitude * latitude)
    return latitude - delta_lat, latitude + delta_lat, longitude - delta_lon, longitude + delta_lon


def distance_square(lat1, lon1, lat2, lon2):
    ''' @var lat1, lon1, lat2, lon2: integer coordinates multiplied by 100000
        @return: square of distance in meters
        delta_lat 1 == delta_dis 1.1m
        delta_lon 1 == delta_dis 1.1m * (1-lat*lat/8100/1e+10)
    '''
    lat_mid = (lat1 + lat2) >> 1
    lat = abs(lat1 - lat2)
    lon = abs(lon1 - lon2) * (81000000000000 - lat_mid * lat_mid) / 81000000000000
    return 121 * (lat * lat + lon * lon) / 100


def distance_haversine(lat1, lon1, lat2, lon2):
    '''
    @var lat1, lon1, lat2, lon2: standar coordinates (NOT multiplied by 100000)
    @return: distance in kilometers
    '''
    R = 6371 #km
    dLat = math.radians(lat2-lat1)
    dLon = math.radians(lon2-lon1)
    lat1 = math.radians(lat1)
    lat2 = math.radians(lat2)
    a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLon/2) * math.sin(dLon/2) * math.cos(lat1) * math.cos(lat2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = R * c
    return d


if __name__ == "__main__":
    log_info("Please!")


